#ifndef MODEL_TURBO
#define MODEL_TURBO


#include <vector>
#include <map>
#include <iostream>
#include <fstream>
#include <sstream>
#include "patch.h"
#include "grid.h"
#include "parallel.h"

//Vector of variables +  overloaded operators
typedef std::vector<double> var_vector;
std::vector<double> operator+=(std::vector<double>& a,const std::vector<double>& b)
{
	if (a.size()!=b.size()) throw 1;
	for (int i = 0; i<a.size(); i++) a[i] += b[i];
	return a;
};

std::vector<double> operator-=(std::vector<double>& a, const std::vector<double>& b)
{
	if (a.size()!=b.size()) throw 1;
	for (int i = 0; i<a.size(); i++) a[i] -= b[i];
	return a;
};

std::vector<double> operator+(const std::vector<double>& a,const std::vector<double>& b)
{		
	var_vector res(a);
	res += b;
	return res;
};

std::vector<double> operator-(const std::vector<double>& a, const std::vector<double>& b)
{		
	var_vector res(a);
	res -= b;
	return res;
};

double operator*(const std::vector<double>& a,const std::vector<double>& b)
{		
	double res = 0;
	for (int i = 0; i<b.size(); i++) res += a[i]*b[i];
	return res;
};

std::vector<double> operator*(const double& a,const std::vector<double>& b)
{		
	var_vector res(b);
	for (int i = 0; i<res.size(); i++) res[i] *= a;	
	return res;
};

std::vector<double> operator*(const std::vector<double>& b,const double& a)
{		
	var_vector res(b);
	for (int i = 0; i<res.size(); i++) res[i] *= a;	
	return res;
};

std::vector<double> operator/(std::vector<double>& b, const double& a)
{		
	var_vector res(b);
	for (int i = 0; i<res.size(); i++) res[i] /= a;	
	return res;
};

//Describes mathematic model of problem
class Model {
public:
	int nv;	//Number of base variables in var_vector
	int nc;	//Number of gas components
	std::vector<double> molarmass; //Molar masses of components
	Grid grid;	//Grid structure
	std::map<int, Cell> intercells;	//Cells that out of grid
	Cell& GetCell(int Global_index); //Accessor to cell data
	Parallel& P;	//Parallel implementation
	std::vector<var_vector> values;	//Array of values in grid cells
	std::map<int, var_vector> outvalues;	//Values acquired from exchanges with other processes (GlobalInd -> value)
	std::vector<Patch*> BCTypes;	//Array of availible boundary condition types
	std::map<int, Patch*> BCCondition;	//Boundary conditions corresponding BCMarker for grid faces
	//Load config and parse parameters	
	int InitBoundary();	//Define all available boundary conditions for the  model
	param_list& ConfigParameters;	//Configuration parameters
	int LoadConfig();		
	//Parallel exchange functions
	int InitParallel();	//Init all parallel operations for the model
	int ExchangeValues();	//Exchange values with other processes
	int GatherValues();	
	var_vector InitConditions(Cell& c);	//Initial conditions for task
	//Accessor for model values
	var_vector GetValue(int Global_index);
	//Constructor for parallel task
	Model(Grid& _grid, Parallel& _P, param_list& lst);
	//Boundary conditions handler
	var_vector GetBoundaryValue(Face& f, var_vector UL);
	var_vector GetBoundaryVal(Face& f, var_vector UL);
	std::vector<Vector> GetBoundaryGradient(Face& f, Cell&c, var_vector U, std::vector<Vector>& grad);
	//Model dependent parameters
	double gamma;	//adiabatic exponent of a gas
	double eps;		//Harten correction parameter
	double CFL;		//CFL value
	double my;		//dynamic Viscosity
	double Pr; 		//Prandtl
	double R;	//Universal gas constant
	//Other functions like jacobians and so on
	//...
	var_vector Fi(var_vector);
	var_vector Gi(var_vector);
	void ToTecPlot(const char*);
	void ToTecPlotHex(const char*);
	void ToTecPlot3Dsurface(const char*, int);
	void ToTecPlot1D(const char*);
	void ToTecPlot2D(const char*, int N, int M);
	void ToTecPlot2(const char*);
	void ToTecPlotSol2DTr(const char *);
	void ToTecPlotSol2DSq(const char *);
};

//Get cell data
Cell& Model::GetCell(int Global_index) {	
	if (P.IsLocal(Global_index)) return grid.cells[grid.cell_index[Global_index]];
	if (intercells.find(Global_index) != intercells.end()) {		
		return intercells[Global_index];
	};	
	throw 1;	
	Cell c;
	return c;
};

//Init model. Get right part of grid and do other preparations
Model::Model(Grid& _g, Parallel& _p, param_list& _lst) : P(_p) , ConfigParameters(_lst) {	
	//Load number of gas components
	std::istringstream ss;
	std::string ncomponents = GetParam(ConfigParameters, "NumberOfComponents")[0];
	ss.str(ncomponents);
	ss>>nc;
	nv = 4 + nc;
	//Init grid and parallel part		
	//Return to each process its own part of grid
	std::vector<bool> IsCell(P.pnodes.size());
	for (int i = 0; i<P.pnodes.size(); i++) IsCell[i] = (P.pnodes[i] == P.rank);					
	grid = _g.SubGrid(IsCell);	
	//Store cells that is on intersection of grids
	intercells.clear();
	for (int i = 0; i<grid.faces.size(); i++) if (!grid.faces[i].isExternal) {
		Face& f = grid.faces[i];
		if (P.IsLocal(f.FaceCell_1) && !P.IsLocal(f.FaceCell_2)) 
			intercells[f.FaceCell_2] = _g.cells[_g.cell_index[f.FaceCell_2]];
		if (!P.IsLocal(f.FaceCell_1) && P.IsLocal(f.FaceCell_2)) 
			intercells[f.FaceCell_1] = _g.cells[_g.cell_index[f.FaceCell_1]];
	};	
	values.resize(grid.cells.size());
};


//Get value int given cell
var_vector Model::GetValue(int Global_index) {
	if (P.IsLocal(Global_index)) return values[grid.cell_index[Global_index]];
	if (outvalues.find(Global_index) != outvalues.end()) {
		//printf("rank %d , value for cell N%d is %lg\n", P.rank, Global_index, outvalues[Global_index][0]);
		return outvalues[Global_index];
	};
	//printf("Hey %d! I wanna cell value %d\n", P.rank, Global_index);
	throw 1;	
	return var_vector(0);
};

//Init parallel interaction (nv must be set)
int Model::InitParallel() {	
	P.AddRequestGroup("values");
	std::vector<std::pair<int, int> > cells_to_send(0);	//what cell to send and where to send	
	std::vector<std::pair<int, int> > cells_to_recv(0); //what cell to recv and from whom to reciev	
	for (int i = 0; i<grid.faces.size(); i++) {
		Face Fi = grid.faces[i];
		if (Fi.isExternal) continue;	//Ignore external faces
		bool isFace1 = P.IsLocal(Fi.FaceCell_1);
		bool isFace2 = P.IsLocal(Fi.FaceCell_2);		
		if ((isFace1) && (!isFace2)) {
			std::pair<int, int> cell_dest;
			cell_dest.first = Fi.FaceCell_1;
			cell_dest.second = P.GetRank(Fi.FaceCell_2);
			cells_to_send.push_back(cell_dest);

			std::pair<int, int> cell_recv;
			cell_recv.first = Fi.FaceCell_2;
			cell_recv.second = P.GetRank(Fi.FaceCell_2);			
			cells_to_recv.push_back(cell_recv);
		};
		if ((!isFace1) && (isFace2)) {
			std::pair<int, int> cell_dest;
			cell_dest.first = Fi.FaceCell_2;
			cell_dest.second = P.GetRank(Fi.FaceCell_1);
			cells_to_send.push_back(cell_dest);			
			
			std::pair<int, int> cell_recv;
			cell_recv.first = Fi.FaceCell_1;
			cell_recv.second = P.GetRank(Fi.FaceCell_1);			
			cells_to_recv.push_back(cell_recv);
		};
	};	
	//Initialize send requests
	for (int i  = 0; i<cells_to_send.size(); i++) {
		//printf("i %ld rank %ld : %ld -> %ld\n", i , P.rank, cells_to_send[i].first, cells_to_send[i].second);
		P.request_groups["values"]->AddSendRequest(cells_to_send[i].first, nv, cells_to_send[i].second, 0);//cells_to_send[i].first);		
	};		
	//Initialize recv requests
	for (int i  = 0; i<cells_to_recv.size(); i++) {
		//printf("i %ld rank %ld : %ld <- %ld\n", i , P.rank, cells_to_recv[i].first, cells_to_recv[i].second);
		P.request_groups["values"]->AddRecvRequest(cells_to_recv[i].first, nv, cells_to_recv[i].second, 0);//cells_to_recv[i].first);		
	};
	//Initialize datastructure for exchange cell values
	outvalues.clear();
	return 0;
};

//Exchange values between processes
int Model::ExchangeValues() {	
	//First we send all we need	
	for (int i = 0; i<P.request_groups["values"]->requests_to_send.size(); i++) {		
		var_vector val = GetValue(P.request_groups["values"]->requests_to_send[i]->cell_index);
		for (int j = 0; j<val.size(); j++) {
			P.request_groups["values"]->requests_to_send[i]->buf[j] = val[j];
		};
	};
	P.request_groups["values"]->SendAll();	
	//Then we recv all needed data
	P.request_groups["values"]->RecvAll();		
	outvalues.clear();
	for (int i = 0; i<P.request_groups["values"]->requests_to_recv.size(); i++){		
		var_vector val(nv);
		for (int j = 0; j<val.size(); j++) {
			val[j] = P.request_groups["values"]->requests_to_recv[i]->buf[j];
		};		
		outvalues[P.request_groups["values"]->requests_to_recv[i]->cell_index] = val;
	};
	//Wait for others to complete
	P.Barrier();
	return 0;
};

//Gather values from all processes
int Model::GatherValues() {
	return 0;
};

//Load all available boundary conditions for the  model
int Model::InitBoundary()
{
	//Here we list all boundary conditions supported by this model
	BCTypes.push_back(new PatchNoSlip());
	BCTypes.push_back(new PatchNoLeak());
	BCTypes.push_back(new PatchSteadyFlow());
	BCTypes.push_back(new PatchRotate());
	BCTypes.push_back(new PatchOutFlowP(gamma));
	BCTypes.push_back(new PatchInFlow(gamma, R));
	BCTypes.push_back(new PatchConstant(gamma, R));
	BCTypes.push_back(new PatchConstantP(gamma,	R));

	//Go through BC params and compose marker -> condition relationship
	std::istringstream ss;
	std::vector<std::string> markers = GetParam(ConfigParameters, "BC");
	for (int i = 0; i<markers.size(); i++) {	
		int bc_marker;
		std::string bc_name;		
		ss.clear(); ss.str(markers[i]); ss >> bc_marker >> bc_name;
		for (int i = 0; i<BCTypes.size(); i++) {
			if (BCTypes[i]->Name == (bc_name)) {
				BCCondition[bc_marker] = BCTypes[i]->Copy();
				BCCondition[bc_marker]->ParseParams(ss);
			};
		};
	};

	return 0;
};

var_vector Model::GetBoundaryValue(Face& f, var_vector UL)
{
	return BCCondition[f.BCMarker]->Apply(f, UL);
};

var_vector Model::GetBoundaryVal(Face& f, var_vector UL)
{
	return BCCondition[f.BCMarker]->ApplyVal(f,UL);
};

std::vector<Vector> Model::GetBoundaryGradient(Face& f, Cell&c, var_vector UL, std::vector<Vector>& grad)
{
	return BCCondition[f.BCMarker]->ApplyVisc(f, c, UL, grad);
};

//Load configuration
int Model::LoadConfig()
{	
	std::istringstream ss;	
	//Get all necessary parameters from list
	std::string tmp = GetParam(ConfigParameters, "MolarMasses")[0];
	ss.clear(); ss.str(tmp);
	while (!ss.eof()) {
		double mm;
		ss >> mm;
		molarmass.push_back(mm);		
	};
	tmp = GetParam(ConfigParameters, "Gamma")[0];
	ss.clear(); ss.str(tmp); ss >> gamma;
	tmp = GetParam(ConfigParameters, "Viscosity")[0];
	ss.clear(); ss.str(tmp); ss >> my;
	tmp = GetParam(ConfigParameters, "Harten")[0];
	ss.clear();ss.str(tmp); ss >> eps;
	tmp = GetParam(ConfigParameters, "Pr")[0];
	ss.clear();ss.str(tmp); ss >> Pr;
	tmp = GetParam(ConfigParameters, "CFL")[0];
	ss.clear();ss.str(tmp); ss >> CFL;
	tmp = GetParam(ConfigParameters, "R")[0];
	ss.clear();ss.str(tmp); ss >> R;
	
	//Acquire initial conditions	
	for (int i = 0; i<values.size(); i++) {
		values[i] = InitConditions(grid.cells[i]);
	};
	//ToTecPlot2DSq("init.dat");
	//Initialize boundary conditions	
	InitBoundary();
	return 0;
};


//Apply init conditions
var_vector Model::InitConditions(Cell& cell)
{
	Vector point = cell.CellCenter;
	std::istringstream ss;	
	std::string InitConditionsType = GetParam(ConfigParameters, "InitConditionType")[0];
	//Recognize particular init conditions type
	if (InitConditionsType == "Simple") {
		if (nc > 1) {
			throw new Exception("Simple init conditions dont fit for multicomponent gas\n");
		};
		std::string tmp = GetParam(ConfigParameters, "InitValues")[0];		
		var_vector val(nv);
		double ro, u , w, v, P, T;
		ss.clear();ss.str(tmp); ss>>u>>w>>v>>P>>T;
		ro = P/(R*T);		
		val[0] = ro;
		val[1] = u*ro;
		val[2] = w*ro;
		val[3] = v*ro;
		val[4] = P/(gamma-1.0) + ro*(u*u+v*v+w*w)/2.0;		
		return val;
	};

	if (InitConditionsType == "Linear") {
		if (nc > 1) {
			throw new Exception("Linear init conditions dont fit for multicomponent gas\n");
		};
		std::string tmp = GetParam(ConfigParameters, "InitValues")[0];		
		var_vector val(nv);
		double ro, u , w, v, P, T;
		ss.clear();ss.str(tmp); ss>>u>>w>>v>>P>>T;
		tmp = GetParam(ConfigParameters, "InitPoint")[0];
		Vector init_point;
		ss.clear();ss.str(tmp); ss>>init_point.x>>init_point.y>>init_point.z;
		tmp = GetParam(ConfigParameters, "GradVx")[0];
		Vector grad_vx;
		ss.clear();ss.str(tmp); ss>>grad_vx.x>>grad_vx.y>>grad_vx.z;
		tmp = GetParam(ConfigParameters, "GradVy")[0];
		Vector grad_vy;
		ss.clear();ss.str(tmp); ss>>grad_vy.x>>grad_vy.y>>grad_vy.z;
		tmp = GetParam(ConfigParameters, "GradVz")[0];
		Vector grad_vz;
		ss.clear();ss.str(tmp); ss>>grad_vz.x>>grad_vz.y>>grad_vz.z;
		tmp = GetParam(ConfigParameters, "GradT")[0];
		Vector grad_T;
		ss.clear();ss.str(tmp); ss>>grad_T.x>>grad_T.y>>grad_T.z;
		tmp = GetParam(ConfigParameters, "GradP")[0];
		Vector grad_P;
		ss.clear();ss.str(tmp); ss>>grad_P.x>>grad_P.y>>grad_P.z;
		u += (point - init_point) * grad_vx;
		w += (point - init_point) * grad_vy;
		v += (point - init_point) * grad_vz;
		T += (point - init_point) * grad_T;
		P += (point - init_point) * grad_P;
		ro = P/(R*T);		
		val[0] = ro;
		val[1] = u*ro;
		val[2] = w*ro;
		val[3] = v*ro;
		val[4] = P/(gamma-1.0) + ro*(u*u+v*v+w*w)/2.0;		
		return val;
	};

	if (InitConditionsType == "Zhukov") {
		if (molarmass.size() != 2) {
			throw new Exception("Zhukov init conditions fit only for 2 component gas\n");
		};
		double Lx = 72e-3;
		double Ly = 150e-3;
		int n_modes;
		std::string tmp = GetParam(ConfigParameters, "InitModes")[0];
		ss.clear();ss.str(tmp); ss>>n_modes;
		tmp = GetParam(ConfigParameters, "InitValues")[0];						
		var_vector val(nv);
		double ro, u ,v , w, P, T;
		double P0;
		double ro1, ro2;
		u = w = v = 0;
		ss.clear();ss.str(tmp); ss>>P0>>ro1>>ro2;
		double y_plane = 100e-3;
		double x0 = 36e-3;
		double y0 = 100e-3;

		//Calculate surface function
		double f = 0;
		const double Pi = 3.1415926535897;
		const int n = 6;
		int m[n] = {0, 3, 3, 2, 4, 2};
		double fi[n];
		int p[n] = {2, 3, 5, 7, 11, 13};
		double k[n];
		for (int i = 0; i<n; i++) {
			fi[i] = 2.0*Pi*m[i]/8;			
			k[i] = (2*Pi/Lx)*p[m[i]];
		};
		
		tmp = GetParam(ConfigParameters, "UniformGravityVector")[0];
		Vector g;
		ss.clear();ss.str(tmp); ss>>g.x>>g.y>>g.z;
		
		double alpha0 = 0.5;
		for (int i = 0; i<n_modes; i++) {
			for (int j = 0; j<n; j++) {
				double x = point.x - x0;
				double y = 0;		
				double a = alpha0 / sqrt(k[i]*k[i]);
				f -= a * cos(k[i]*x + fi[i]);
				break;
			};
		};
						
		double B;
		double c2;
		double dy = (point.y - y0); 
		double h =  dy + f;
		v = 0 * dy;
		if ( dy > 0 ) {
			//P = P0 + g.y * ro1 * dy / 2;
		} else {
			//P = P0 + g.y * ro2 * dy / 2;	
		};
		if ( h > 0) {			
			ro = ro1;
			c2 = gamma*P0/ro;
			B = gamma * g.y / c2;
			P = P0 * exp( B * dy );
			ro *= exp( B * dy );
			//P = P0 + ro*g.y*h;
			//ro *= P/P0;
		} else {
			ro = ro2;
			c2 = gamma*P0/ro;
			B = gamma * g.y / c2;
			P = P0 * exp( B * dy );
			ro *= exp( B * dy );
			//P = P0 + ro * g.y * h;
			//ro *= P/P0;
		};		
		val[0] = ro;
		val[1] = u*ro;
		val[2] = v*ro;
		val[3] = w*ro;
		val[4] = P/((gamma-1.0)) + ro*(u*u+v*v+w*w)/2.0;		
		return val;
	};
	var_vector val(nv);
	return val;
};

//Output solution to tecplot
void Model::ToTecPlotHex(const char *fname)
{	
	std::ofstream ofs(fname);			
	//TO DO unify
	ofs<<"VARIABLES= \"X\", \"Y\", \"Z\", \"Rho\", \"u\", \"v\", \"w\", \"T\", \"P\", \"pnode\"\n";
	ofs<<"ZONE T=\"D\"\n";
	ofs<<"N=" << grid.nodes.size() << ", E=" << grid.cells.size() <<", F=FEBLOCK, ET=BRICK\n";	
	ofs<<"VARLOCATION = (NODAL, NODAL, NODAL, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED)\n";
	//X
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs<<grid.nodes[i].P.x<<"\n";
	};
	//Y
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs<<grid.nodes[i].P.y<<"\n";
	};
	//Z
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs<<grid.nodes[i].P.z<<"\n";
	};
	
	//Rho
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][0]<<"\n";
	};
	//u
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][1]/values[i][0]<<"\n";
	};
	//v
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][2]/values[i][0]<<"\n";
	};
	//w
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][3]/values[i][0]<<"\n";
	};
	
	//T
	for (int i = 0; i<values.size(); i++) {
		double ro = values[i][0];
		double vx = values[i][1]/values[i][0];
		double vy = values[i][2]/values[i][0];
		double vz = values[i][3]/values[i][0];
		double E = values[i][4]/values[i][0];
		double T = (gamma-1.0) * (E - (vx*vx+vy*vy+vz*vz)/2.0) / R;
		ofs<<T<<"\n";
	};
	
	//P
	for (int i = 0; i<values.size(); i++) {
		double ro = values[i][0];
		double vx = values[i][1]/values[i][0];
		double vy = values[i][2]/values[i][0];
		double vz = values[i][3]/values[i][0];
		double E = values[i][4]/values[i][0];
		double P = (gamma-1.0) * ro * (E - (vx*vx+vy*vy+vz*vz)/2.0);
		ofs<<P<<"\n";
	};
	
	//pnode
	for (int i = 0; i<values.size(); i++) {
		ofs<<P.pnodes[i]<<"\n";
	};
	//Connectivity list
	for (int i = 0; i<grid.cells.size(); i++) {	
		ofs<<grid.node_index[grid.cells[i].Nodes[0]]+1<<" "
			<<grid.node_index[grid.cells[i].Nodes[1]]+1<<" "
			<<grid.node_index[grid.cells[i].Nodes[2]]+1<<" "
			<<grid.node_index[grid.cells[i].Nodes[3]]+1<<" "
			<<grid.node_index[grid.cells[i].Nodes[4]]+1<<" "
			<<grid.node_index[grid.cells[i].Nodes[5]]+1<<" "
			<<grid.node_index[grid.cells[i].Nodes[6]]+1<<" "
			<<grid.node_index[grid.cells[i].Nodes[7]]+1<<"\n";
	};
	ofs.close();	
	return;
};
//Output 3D-surface distribution of values v for quadrilateral mesh
void Model::ToTecPlot3Dsurface(const char *fname, int BM)
{
	int k=0;
	std::vector<int> BN_ind(0);
	std::map<int, int> B_node_map;	//Map from Global_index of node to index in array BN_ind
	for(int i=0; i<grid.nodes.size(); i++)
	{
		if(grid.nodes[i].BMarker==BM)
		{
			BN_ind.push_back(i);
			B_node_map[grid.nodes[i].GlobalIndex] = k;
			k++;
		};
	};
	std::vector<int> BF_ind(0);
	for(int i=0; i<grid.faces.size(); i++)
	{
		if((grid.faces[i].isExternal==1)&&(grid.faces[i].BCMarker==BM))
		{
			BF_ind.push_back(i);
		};
	};
	std::ofstream ofs(fname);			
	//TO DO unify
	ofs<<"VARIABLES= \"X\", \"Y\", \"Z\", \"T\"\n";
	ofs<<"ZONE T=\"D\"\n";
	ofs<<"N=" << BN_ind.size() << ", E=" << BF_ind.size() <<", F=FEBLOCK, ET=QUADRILATERAL\n";	
	ofs<<"VARLOCATION = (NODAL, NODAL, NODAL, CELLCENTERED)\n";
	//X
	for (int i = 0; i<BN_ind.size(); i++) {
		ofs<<grid.nodes[BN_ind[i]].P.x<<"\n";
	};
	//Y
	for (int i = 0; i<BN_ind.size(); i++) {
		ofs<<grid.nodes[BN_ind[i]].P.y<<"\n";
	};
	//Z
	for (int i = 0; i<BN_ind.size(); i++) {
		ofs<<grid.nodes[BN_ind[i]].P.z<<"\n";
	};
	//T
	for (int i = 0; i<BF_ind.size(); i++) {
		Face B_face = grid.faces[BF_ind[i]];
		int val_ind = grid.cell_index[B_face.FaceCell_1];
		double ro = values[val_ind][0];
		double vx = values[val_ind][1]/ro;
		double vy = values[val_ind][2]/ro;
		double vz = values[val_ind][3]/ro;
		double E = values[val_ind][4]/ro;		//full energy divided by mass
		double T = (gamma-1.0) * (E - (vx*vx+vy*vy+vz*vz)/2.0) / R;	//R is specific
		ofs<<T<<"\n";
	};
	//Connectivity list
	for (int i = 0; i<BF_ind.size(); i++) {
		Face our_face = grid.faces[BF_ind[i]];
		ofs << B_node_map[our_face.FaceNodes[0]]+1 << " "
			<< B_node_map[our_face.FaceNodes[1]]+1 << " "
			<< B_node_map[our_face.FaceNodes[2]]+1 << " ";
			if(our_face.FaceNodes.size()==4)
		ofs	<< B_node_map[our_face.FaceNodes[3]]+1 << "\n";
			else	//if some faces are triangels
		ofs	<< B_node_map[our_face.FaceNodes[2]]+1 << "\n";
	};
	ofs.close();	
	return;
};

void Model::ToTecPlot1D(const char *fname) {
	std::ofstream ofs(fname);
	//TO DO unify
	ofs<<"VARIABLES= \"X\", \"Rho\", \"P\", \"U\", \"v\", \"w\", \"e\"\n";
	ofs<<"ZONE T=\"D\"\n";
	ofs<<"I=" << grid.cells.size() <<"\n";
	for (int i = 0; i<grid.cells.size(); i++) {
		ofs << grid.cells[i].CellCenter.x << " ";
		ofs << values[i][0] << " ";
		ofs << (gamma-1.0)*(values[i][4] -values[i][1]*values[i][1]/(2.0*values[i][0])) << " "; //P
		ofs << values[i][1]/values[i][0] << " ";
		ofs << values[i][2]/values[i][0] << " ";
		ofs << values[i][3]/values[i][0] << " ";
		ofs << values[i][4]/values[i][0] << "\n";
	};
	ofs.close();
	return;
};

//For regular mesh
void Model::ToTecPlot2D(const char *fname, int N, int M) {
	std::ofstream ofs(fname);
	//TO DO unify
	ofs<<"VARIABLES= \"X\", \"Y\", \"Rho\", \"P\", \"Rho*U\", \"Rho*v\", \"Rho*w\", \"Rho*e\"\n";
	ofs<<"ZONE T=\"D\"\n";
	ofs<<"I=" << N-1 <<" J=" << M-1 <<" F=POINT\n";
	for (int i = 0; i<N-1; i++) {
		for (int j = 0; j<M-1; j++) {
			int ind = i +j*(N-1);
			ofs << grid.cells[ind].CellCenter.x << " ";
			ofs << grid.cells[ind].CellCenter.y << " ";
			ofs << values[ind][0] << " ";
			ofs << (gamma-1.0)*(values[ind][4] -values[ind][1]*values[ind][1]/(2.0*values[ind][0])) << " "; //P
			ofs << values[ind][1] << " ";
			ofs << values[ind][2] << " ";
			ofs << values[ind][3] << " ";
			ofs << values[ind][4] << "\n";
		};
	};	
	ofs.close();
	return;
};

//Output solution to tecplot
void Model::ToTecPlot(const char *fname)
{
	//Debug
	std::vector<int> badcells(0);
	for (int i = 0; i<grid.cells.size(); i++) {
		if (fabs(values[i][0] - 1.0) > 1e-8) {
			badcells.push_back(i);
		};
	};
	std::ofstream ofs2("bad.dat");
	std::ofstream ofs(fname);
	//TO DO unify
	ofs<<"VARIABLES= \"X\", \"Y\", \"Z\", \"Rho\", \"Rho*U\", \"Rho*v\", \"Rho*w\", \"Rho*e\"\n";
	ofs<<"ZONE T=\"D\"\n";
	ofs<<"N=" << grid.nodes.size() << ", E=" << grid.cells.size() <<", F=FEBLOCK, ET=TETRAHEDRON\n";	
	ofs<<"VARLOCATION = (NODAL, NODAL, NODAL, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED)\n";
	//X
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs<<grid.nodes[i].P.x<<"\n";
	};
	//Y
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs<<grid.nodes[i].P.y<<"\n";
	};
	//Z
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs<<grid.nodes[i].P.z<<"\n";
	};
	
	//Rho
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][0]<<"\n";
	};
	//Rho*u
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][1]<<"\n";
	};
	//Rho*v
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][2]<<"\n";
	};
	//Rho*w
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][3]<<"\n";
	};
	//Rho*e
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][4]<<"\n";
	};
	
	//Connectivity list
	for (int i = 0; i<grid.cells.size(); i++) {	
		ofs<<grid.cells[i].Nodes[0]+1<<" "
			<<grid.cells[i].Nodes[1]+1<<" "
			<<grid.cells[i].Nodes[2]+1<<" "
			<<grid.cells[i].Nodes[3]+1<<"\n";
	};
	ofs.close();
	//Debug
	ofs2<<"VARIABLES= \"X\", \"Y\", \"Z\", \"Rho\", \"Rho*U\", \"Rho*v\", \"Rho*w\", \"Rho*e\"\n";
	ofs2<<"ZONE T=\"D\"\n";
	ofs2<<"N=" << grid.nodes.size() << ", E=" << badcells.size() <<", F=FEBLOCK, ET=TETRAHEDRON\n";
	ofs2<<"VARLOCATION = (NODAL, NODAL, NODAL, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED)\n";
	//X
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs2<<grid.nodes[i].P.x<<"\n";
	};
	//Y
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs2<<grid.nodes[i].P.y<<"\n";
	};
	//Z
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs2<<grid.nodes[i].P.z<<"\n";
	};
	
	//Connectivity list	
	//Rho
	for (int i = 0; i<badcells.size(); i++) {
		ofs2<<values[badcells[i]][0]<<"\n";
	};
	//Rho*u
	for (int i = 0; i<badcells.size(); i++) {
		ofs2<<values[badcells[i]][1]<<"\n";
	};
	//Rho*v
	for (int i = 0; i<badcells.size(); i++) {
		ofs2<<values[badcells[i]][2]<<"\n";
	};
	//Rho*w
	for (int i = 0; i<badcells.size(); i++) {
		ofs2<<values[badcells[i]][3]<<"\n";
	};
	//Rho*e
	for (int i = 0; i<badcells.size(); i++) {
		ofs2<<values[badcells[i]][4]<<"\n";
	};
	for (int i = 0; i<badcells.size(); i++) {
		ofs2<<grid.cells[badcells[i]].Nodes[0]+1<<" "
			<<grid.cells[badcells[i]].Nodes[1]+1<<" "
			<<grid.cells[badcells[i]].Nodes[2]+1<<" "
			<<grid.cells[badcells[i]].Nodes[3]+1<<"\n";
	};
	ofs2.close();
	return;
};


//Second output
void Model::ToTecPlot2(const char *fname)
{
	std::ofstream ofs(fname);
	//TO DO unify
	ofs<<"VARIABLES= \"X\", \"Y\", \"Z\", \"Rho\", \"Rho*U\", \"Rho*v\", \"Rho*w\", \"Rho*e\"\n";
	ofs<<"ZONE T=\"D\"\n";	
	
	//Rho
	for (int i = 0; i<values.size(); i++) {
		ofs<<grid.cells[i].CellCenter.x<<" "
			<<grid.cells[i].CellCenter.y<<" "
			<<grid.cells[i].CellCenter.z<<" "
			<<values[i][0]<<" "
			<<values[i][1]<<" "
			<<values[i][2]<<" "
			<<values[i][3]<<" "			
			<<values[i][4]<<"\n";
	};
	
	ofs.close();
	return;
};



void Model::ToTecPlotSol2DTr(const char *fname_sol)
{
	std::ofstream ofs(fname_sol);
	//TO DO unify
	ofs<<"VARIABLES= \"X\", \"Y\", \"Rho\", \"U\", \"v\", \"w\", \"e\"\n";
	ofs<<"ZONE T=\"D\"\n";
	ofs<<"N=" << grid.nodes.size() << ", E=" << grid.cells.size() <<", F=FEBLOCK, ET=TRIANGLE\n";	
	ofs<<"VARLOCATION = (NODAL, NODAL, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED)\n";
	//X
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs<<grid.nodes[i].P.x<<"\n";
	};
	//Y
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs<<grid.nodes[i].P.y<<"\n";
	};
		
	//Rho
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][0]<<"\n";
	};
	//u
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][1]/values[i][0]<<"\n";
	};
	//v
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][2]/values[i][0]<<"\n";
	};
	//w
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][3]/values[i][0]<<"\n";
	};
	//e
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][4]/values[i][0]<<"\n";
	};
	
	//Connectivity list
	for (int i = 0; i<grid.cells.size(); i++) {	
		ofs<<grid.cells[i].Nodes[0]<<" "
			<<grid.cells[i].Nodes[1]<<" "
			<<grid.cells[i].Nodes[2]<<"\n";
	};
	ofs.close();
	return;
};

//Output solution to tecplot	
void Model::ToTecPlotSol2DSq(const char *fname_sol)
{
	std::ofstream ofs(fname_sol);
	//TO DO unify		
	ofs<<"VARIABLES= \"X\", \"Y\", \"Rho\", \"U\", \"v\", \"w\", \"e\", \"T\", \"P\", \"proc\"\n";
	ofs<<"ZONE T=\"D\"\n";
	ofs<<"N=" << grid.nodes.size() << ", E=" << grid.cells.size() <<", F=FEBLOCK, ET=QUADRILATERAL\n";
	ofs<<"VARLOCATION = (NODAL, NODAL, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED)\n";
	//X
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs<<grid.nodes[i].P.x<<"\n";
	};
	//Y
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs<<grid.nodes[i].P.y<<"\n";
	};

	//Rho
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][0]<<"\n";
	};
	//u
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][1]/values[i][0]<<"\n";
	};
	//v
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][2]/values[i][0]<<"\n";		
	};
	//w
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][3]/values[i][0]<<"\n";
	};
	//e
	for (int i = 0; i<values.size(); i++) {
		ofs<<values[i][4]/values[i][0]<<"\n";
	};
	//T
	for (int i = 0; i<values.size(); i++) {
		double ro = values[i][0];
		double vx = values[i][1]/values[i][0];
		double vy = values[i][2]/values[i][0];
		double vz = values[i][3]/values[i][0];
		double E = values[i][4]/values[i][0];
		double T = (gamma-1.0) * (E - (vx*vx+vy*vy+vz*vz)/2.0) / R;
		ofs<<T<<"\n";
	};
	//P
	for (int i = 0; i<values.size(); i++) {
		double ro = values[i][0];
		double vx = values[i][1]/values[i][0];
		double vy = values[i][2]/values[i][0];
		double vz = values[i][3]/values[i][0];
		double E = values[i][4]/values[i][0];
		double P = (gamma-1.0) * ro * (E - (vx*vx+vy*vy+vz*vz)/2.0);
		ofs<<P<<"\n";
	};
	//TO DO	
	for (int i = 0; i<grid.cells.size(); i++) {		
		ofs<<P.pnodes[i]<<"\n";
	};	
	//automodel
	std::ofstream test("test.dat");

	for (int i = 0; i<values.size(); i++) {
		double x = grid.cells[i].CellCenter.x;
		double y = grid.cells[i].CellCenter.y;
		double res = 0;
		if ((x>=1.25)&&(x<1.26)) {
			res = y*sqrt((1.0)/(0.002*(x - 0.25)));
			test<<res<<" "<<values[i][1]<<"\n";
		};
	};

	//Connectivity list
	for (int i = 0; i<grid.cells.size(); i++) {
		ofs<<grid.cells[i].Nodes[0]<<" "
			<<grid.cells[i].Nodes[1]<<" "
			<<grid.cells[i].Nodes[2]<<" "
			<<grid.cells[i].Nodes[3]<<"\n";
	};
	test.close();
	ofs.close();
	return;
}

#endif
